% 周期矩形脉冲信号的傅里叶变换与逆变换（矢量法与矩阵法效率对比）
clear; clc; close all;

% 参数设置
T0 = 2;         % 信号周期
tau = 1;        % 脉冲宽度
N_period = 5;   % 展示周期数
fs = 1000;      % 采样频率
N = T0*fs*N_period; % 总采样点数

t = linspace(-N_period/2*T0, N_period/2*T0, N); % 时间轴
x = zeros(size(t));
for k = -ceil(N_period/2):ceil(N_period/2)
    x = x + double(abs(t - k*T0) <= tau/2);
end

% 绘制周期矩形脉冲信号
figure;
plot(t, x, 'LineWidth', 2);
xlabel('时间 t'); ylabel('幅值');
title('周期矩形脉冲信号');
grid on;
% saveas(gcf, './figure/signal_exp2_fig1.png');

% 傅里叶变换参数
L = length(x);
f = (-L/2:L/2-1)*(fs/L); % 频率轴

%% 1. 矢量法实现DFT和IDFT
fprintf('--- 矢量法 ---\n');
tic;
X_vec = fftshift(fft(x));
t_vec = toc;
tic;
x_rec_vec = ifft(ifftshift(X_vec));
t_vec_inv = toc;

% 打印矢量法傅里叶变换前10个复数结果
fprintf('矢量法傅里叶变换结果（前10个复数）：\n');
disp(X_vec(1:10));

figure;
subplot(2,1,1);
plot(f, abs(X_vec), 'LineWidth', 1.5);
xlabel('频率 f'); ylabel('|X(f)|'); title('矢量法-幅度谱'); grid on;
subplot(2,1,2);
plot(t, real(x_rec_vec), 'r'); hold on; plot(t, x, 'k--');
legend('重构信号','原信号'); xlabel('时间 t'); ylabel('幅值'); title('矢量法-逆变换重构'); grid on;
% saveas(gcf, './figure/signal_exp2_fig2.png');

%% 2. 矩阵法实现DFT和IDFT
fprintf('--- 矩阵法 ---\n');
tic;
n = 0:L-1;
k = n';
WN = exp(-1j*2*pi/L);
F = WN.^(k*n); % DFT矩阵
X_mat = (F*x.').';
t_mat = toc;
tic;
F_inv = exp(1j*2*pi/L).^(k*n)/L; % IDFT矩阵
x_rec_mat = (F_inv*X_mat.').';
t_mat_inv = toc;

% 打印矩阵法傅里叶变换前10个复数结果
fprintf('矩阵法傅里叶变换结果（前10个复数）：\n');
disp(X_mat(1:10));

figure;
subplot(2,1,1);
plot(f, abs(fftshift(X_mat)), 'LineWidth', 1.5);
xlabel('频率 f'); ylabel('|X(f)|'); title('矩阵法-幅度谱'); grid on;
subplot(2,1,2);
plot(t, real(x_rec_mat), 'r'); hold on; plot(t, x, 'k--');
legend('重构信号','原信号'); xlabel('时间 t'); ylabel('幅值'); title('矩阵法-逆变换重构'); grid on;
% saveas(gcf, './figure/signal_exp2_fig3.png');

%% 3. 打印效率对比
fprintf('矢量法DFT用时: %.6f 秒\n', t_vec);
fprintf('矢量法IDFT用时: %.6f 秒\n', t_vec_inv);
fprintf('矩阵法DFT用时: %.6f 秒\n', t_mat);
fprintf('矩阵法IDFT用时: %.6f 秒\n', t_mat_inv);
